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Abstract 

We compare the Rubinstein-Duke model for reptation to a model where the boundary dynamics 
is modified by calculating the viscosity of polymer melts. The question is investigated whether 
the viscosity is determined by details of the dynamics of the polymer ends or by the stretching of 
the polymer. To this end the dependence of the viscosity on the particle density of the lattice gas 
models which can be identified by the stretching is determined. We show that the influence of the 
stretching of the polymer on the absolute value of the viscosity in the scaling limit of of very long 
chains is much bigger than the influence of the boundary dynamics, whereas the corrections of the 
scaling of the viscosity depends significantly on the details of the boundary dynamics. 
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I. INTRODUCTION 



A persistent problem in polymer science is the theoretical description of the mass de- 
pendence of the viscosity of polymer melts of linear chains.— A model which recently lead 
to some progress^jM is the Rubinstein-Duke (RD) model,i^ which is a discrete model that 
maps the reptation dynamics of a three dimensional polymer onto a lattice-gas. This model 
takes into account tube length fluctuations as well as the effect of a possible external electric 
field which acts on a charged polymer. 

In the reptation model the stretching of the polymer is caused by an entropic tensile force 
acting on the chain ends.— In the Rubinstein model this force is parameterized by the lattice 
dimension which provides the rates for adding and removing new segments to the tube. 
Adding (resp. removing) new segments is modeled in the lattice-gas picture by the cre- 
ation (resp. annihilation) of particles at the boundaries. The stretching of the polymer 
corresponds to the particle density in the bulk. Although the real spatial dimension of the 
polymer network (naturally d = 3) determines the possibilities of the polymer to stretch 
and by this it is directly related to the creation rate, this rate is not only influenced by the 
dimension because the dynamics of the boundary segments of a polymer is affected by more 
mechanisms than pure geometry. 

Rubinstein investigated the influence of the creation and annihilation rates on the viscosity 
which is proportional to the relaxation time and showed that the magnitude of the viscosity 
increases with increasing creation rates. Due to particle number conservation in the bulk 
the particle density is provided by the boundary dynamics - none the less the density is a 
property of the bulk and different boundary dynamics may lead to the same density. This 
feature of the model corresponds to a realistic property of the polymer: If the behavior of the 
boundary segments is changed, the stretching of the polymer will not necessarily be affected. 
In order to investigate whether the boundary dynamics or the particle density determines 
the relaxation time we consider two models: The standard Rubinstein-Duke model and a 
modified version. The model is only modified at the boundaries in a way that in the bulk 
the particle density and thus the entropic force remains unchanged.— 

We will show that the particle density is the crucial property, so that it is reasonable to 
specify the density dependence of the relaxation time explicitly. As a hypothesis we convert 
the known density dependence of the diffusion constant to an expression for the relaxation 
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time. This hypothesis is confirmed by our numerical data implying that we have on the one 
hand determined the density dependence and on the other hand verified the conversion of 
the diffusion constant to the relaxation time. However, we also show that microscopic details 
of chain end dynamics have an appreciable impact on the absolute value of the viscosity and 
its scaling. 

II. MODEL 

A. Standard Rubinstein Duke model 

Rubinstein assumes that the constraints of the other polymers divide space into cells 
which form a (i-dimensional regular cubic lattice. The polymer occupies a series of adjacent 
cells, the "primitive path". It is not possible for the polymer to traverse the edges of the 
cells (in two dimensions: the lattice points) so that only the ends of the polymer can enter 
new cells. The polymer is divided into segments whose lengths are of the order of the 
lattice constant, the number of segments is proportional to the length of the polymer or 
the molecular weight. The ends of the segments are labeled with imaginary particles called 
"reptons". The orientation of the lattice is introduced by Duke in a way that the electric 
field is diagonal to the lattice, i.e. in three dimensions the (lll)-direction, which is relevant 
e.g. for the investigation of electrophoresis. In this paper we will only consider the limit of 
zero field - nevertheless we will keep a reference axis along which the displacement of the 
polymer chain whole can be monitored.^*"^ 

A repton is allowed to jump into an adjacent cell according to the following rules: 

1. The reptons in the bulk are only allowed to jump along the primitive path. 

2. No cell in the interior of the primitive path may be left empty. 

3. The ends move freely provided rule 2 is respected. If an end repton occupies the cell 
alone, it can only retract in the cell of the adjacent repton. If the adjacent repton is 
in the same cell, the end repton may enter any of the 2d surrounding cells. Reptons 
in the bulk jump with the same probability as reptons at the end into occupied cells. 



3 




FIG. 1: a) The repton model in two dimensions: the circles represent the reptons; the primitive 
path is marked by the bold lines, b) Projection onto one dimension. 

Rule 1 refers to the diffusion of defects and ensures that the polymer does not traverse the 
topological constraints which are represented by the edges of the cells in the model. Rule 
2 refers to the connectivity of the polymer and is motivated by the fact that the segments 
are of the size of the lattice constant and are not allowed to stretch. Finally rule 3 reflects 
the fact that there are more free adjacent cells for an end repton than occupied ones - this 
introduces the entropic tensile force. 

As the shape of the primitive path is not affected by the movement of the polymer and only 
the ends are created or annihilated this model can be mapped onto one dimension. To this 
end a particle representation is used: Starting on one end, a segment into the direction of 
the reference axis is identified with an A-particle, against the axis with a 5-particle and 
finally a segment within one cell by a vacancy (Fig. 

These particles are residing on a chain whose number of sites L is the number of segments, 
thus the number of reptons minus one. As for each site three states {A,B,0} are possible 
this model can be identified with a quantum spin-one chain . 

The one-dimensional dynamics of the particles is the following: In the bulk. A- and B- 
particles interchange freely with vacancies 0, but may not traverse each other. By these 
rules it is guaranteed that the shape of the primitive path is conserved. 
At the boundaries, particles may be annihilated with rate (3, i.e. transformed into vacan- 
cies, or created with rate a, i.e. the transformation of vacancies into particles. Within the 
approach of Rubinstein the annihilation is a process with the same rate as the hopping in 
the bulk which is a basic process (rate is assigned to 1), i.e. (3 = 1. The creation of particles 
is (i-times more probable, because at the boundaries there are 2d cells to enter, d along and 
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d against the direction of the field, while in the bulk there is only one cell to enter. Hence 
according to Rubinstein, a = d. 

In this way a stochastic interacting particle system on a one dimensional chain has been 
defined. As the transitions are independent of the previous history this process is Marko- 
vian. A convenient way to describe the process mathematically is the quantum Hamiltonian 
formalism which we will present here briefiy, for details see Ref. [ill . 

Assigning to the particles and the vacancy one of the three-dimensional basis vectors each, a 
configuration of the spin chain can be represented by the tensor product \ri) of these vectors. 
The probability to be in state \ri) at time t is labeled by -Pry(t). These probabilities of the 
individual states can be combined to a vector: \P{t)) = ^ Pri{t) \v)- Due to the conservation 
of probability the entries of the vector \P{t)) sum up to 1 at any time t. With this definition 
the master equation can be written as: 

j^\P{t)) = -H\P{t)), 

with the stochastic generator 

Here Wr^'^n is the transition probability from state \v() to |?7). In other words the off diagonal 
elements of the matrix H. are the negative transition rates between the respective states and 
the diagonal elements are the sum of the rates leading away from the respective state. 
The creation operators and 6^ are defined by |0) = \A) and iS^ |0) = acting on 
any other state yields zero. The annihilation operators a and h are defined by a \A) = \0) 
and b \B) = |0), again acting on any other state yields zero. Finally we define the number 
operators ra^ = a^a, = h'^h and = 1 — n"^ — n^. By these definitions the stochastic 
generator H of the RD model reads!^ 

L-l 

i7 = 6i(rf,l) + 6L(rf,l) + ^u„ (1) 

n=l 

with 
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As mentioned in the introduction we are interested in the influence of the stretching of the 
polymer on the viscosity of the polymer melt. To this end we present how to determine the 
stretching and then how to calculate the viscosity. 

Each particle contributes a length of the order of the lattice constant c to the primitive path, 
so that the length of the primitive path equals Nc, where is the total number of particles 
on the chain. The stretching of the polymer corresponds to the length to the primitive path 
divided by the number of bonds L - or in the particle picture to the particle density p. The 
particle density p is determined by the creation rate a and the annihilation rate (3. In the 
bulk the density on one site changes only due to the hopping process: 

^(rifc) = ((1 - nk){nk+i + Uk-i)) 

- {{nk){l - rik+i + I - Uk-i)) (2) 
= {rik+i) + {uk-i) - 2{nk). 

Here, k = 2, . . . , L — 1 and Uk = + ■ The first term on the right hand side equals 
the gain of particle which is proportional to the occupation of the site with vacancies and 
proportional to the occupation of the adjacent sites with particles. Accordingly, the second 
term describes the loss of particles. Supplementary at the boundaries k = {1,L} we have 
to take creation and annihilation into account: 

^{uk) = ((1 - nk)ink±i)) - ((nfc)(l - nfcii)) 
at 

+ 2a{l-nk) - P{nk) (3) 
= K±i) - (nk) + 2a- {2a + (3){nk). 

Here, k G {1, L} and the index k + 1 for k = 1 respectively k — 1 for k = L. We have to 
consider the rate 2a because both A- and 5-particles are counted by nk- As we are only 
interested in the stationary state we presume d/dt (nl) = for A; = 1, . . . , L. Eq. ^ then 
yields 

(nl) = a + bk, 

with constants a and b. Finally, Eq. (jS)) determines a = 2a /{2a + f3) and 6 = 0. 
So the stationary total particle number A^* = '^k=i{''^k) calculated to^^ 



and for the particle density p* = N*/L we get 

In the RD model one therefore has p* = -^^^ The viscosity is proportional to the longest 
relaxation time of the stochastic generator^ H which is the inverse energy gap of the spin 
chain. A very efficient algorithm to calculate the lowest excitations of quantum spin chains 
is the density matrix renormalization group (DMRG) algorithm .^^i^^ For the standard RD 
model the usefulness of this algorithm has been demonstratedj-^*^ so that we employ this 
algorithm as well for the modifications of the RD model. The DMRG algorithm is a nu- 
merical method to diagonalize huge matrices whose systematic is known. This is done by 
reducing the Hilbert space to the most important states which are chosen by calculating 
their probability using the density matrix. 



B. Reservoir at the boundaries 



In order to investigate the infiuence of the boundary dynamics the following modification 
of the model is considered: In the bulk the dynamics of the particles remains unchanged. At 
the boundaries the creation and annihilation rates are chosen as if the chain was in contact 
with a particle reservoir, as described in Ref. 0. In the reservoir the density of A-particles 
is labeled with p/2, the density of 5-particles is the same, so that the total density is p. 
The creation of particles at the boundary is identified with the hopping of a particle from 
a boundary site of the reservoir onto a boundary site on the chain. As the probability to 
find an A- or 5-particle on a site of the reservoir is p/2 each, the rate for the creation 
process has to be chosen as p/2. Accordingly, the annihilation of a particle corresponds to 
the hopping of a particle from the boundary of the chain into the reservoir. This hopping 
is only possible when the boundary site of the reservoir is empty which happens with the 
probability 1 — p. Therefore, this is the annihilation rate. Using Eq. (ji} it is easily calculated 
that this boundary dynamics provides a particle density p in the bulk. 

If we want to provide the same particle density as in the RD model with lattice dimension 
d, we therefore have to choose the creation rate a' = d/{2d+ 1) and annihilation rate 
(3' = l/{2d+ 1). The modified rates differ from the previous ones by a factor l/{2d+ 1). 
A physical interpretation of this modification of the model is that the dynamics of the 
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FIG. 2: Viscosity in dependence of the chain length - comparison of the different boundary dy- 
namics. 

boundary reptons is slowed down, as the rates for creation and annihilation are reduced. 
Experimentally this could be seen for polymers with bulky end groups. 

III. CALCULATIONS AND RESULTS 
A. Comparison of boundary dynamics 

We calculate in the standard RD model the relaxation times for d = 3, d = 1, d = 0.5, 
and d = 0.3. Small values of d seem to be rather artificial - however, as the entanglement 
network is not regular cubic but instead contains also very small pores, one may investigate 
smaller effective dimensions than d = 3. An illustrative argument is that it is much easier 
to pull a thread out of the eye of a needle than to pass it through — by this the retracting 
of a polymer should be more probable than entering a new cell. 

For d = 3 one has a density p = 6/7 so that the length of the primitive path is ld=3 = 0.86cL 
with the lattice constant c and number of bonds L - so the polymer is nearly stretched. For 
the other dimensions the length of the primitive path calculates to ld=i = 0.67cL (p = 2/3), 
ld=o,5 = 0.5cL (p = 1/2), and ld=o.3 = 0.38cL (p = 3/8), the length decreases clearly. 
To compare the standard RD model with the modified boundary dynamics, we compute 
the relaxation times in the modified model with the same particle densities. The results for 
both models are plotted in Fig. |21 
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With decreasing d the relaxation time decreases in both models. This can be explained 
by the fact that a shorter tube relaxes faster. In the particle picture this is justified by 
considering that the relaxation process is the traversion of a particle from one end of the 
chain to the other. As for a lower particle density the diffusion is faster this process takes 
less time. 

One can observe two qualitative results: Firstly compared to the RD model the relaxation 
times of the reservoir boundary dynamics are shifted to higher values. Secondly it is obvious 
that the influence of the particle density is much stronger than the influence of the boundary 
dynamics. 

B. Density dependence of the relaxation time 

After having verified in the last section that the particle density is crucial for the relax- 
ation time, in this section the dependence is quantified. To this end a method for computing 
the diffusion constant from the relaxation time is presented. As the density dependence of 
the diffusion constant is knowoi^ii^ this leads to a hypothesis for the density dependence of 
the relaxation time. By verifying this hypothesis not only the density dependence is found 
but also the validity of the computation of the diffusion constant is proved. 
As the RD model is a projection to a specific axis, the diffusion constant along this axis will 
be calculated which differs only by a constant factor from the three dimensional diffusion 
constant. The projected end to end length of the polymer is the difference of the number 
of A- and 5-particles, \Na — Nb\- The one dimensional diffusion can be identified with a 
random walk: After the relaxation time r the polymer has moved by a distance \Na — Nb\ 
in or against the direction of the field, i.e. 



with a„ = ± \Na — Nb\- Assuming that the random walker starts at x(0) = the following 
relation holds (for sufficiently large times t so that the number of steps is well described by 



x{t + r) 



x{t) + a. 



K = t/T): 



K=t/T 




n=l 
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which yields 



K = t/T 



n=l 



because the a„ are uncorrelated: {anCim) = ^n/mC^n- Thus calculating the diffusion constant 
yields 

D 



{{Na-NbY) 



2r 



As described in Ref. |11| the expectation value is calculated by 



{{Na - Nef) = {s\ {Na - Nsf 



with the "summation vector" (s| = (1,1,...) and the stationary state \P*)- With a prod- 
uct ansatz the stationary state can be computed tc^- (here with generalized creation and 
annihilation rates a respectively 13) 



\P*) 



/5 



2a + (3 



so that 



{{Na - Nb?) = L- 



13 
1 

a 

2a 



2a + (3 

Thus we get for the diffusion constant the expression 

Lp 

D = —. 

2t 



Lp. 



(5) 



In Ref. 



. HorQit i 



is shown that the density dependence of the diffusion constant is 



(6) 



The term K{L) corrects finite size effects of the length dependence. Thus the density 
dependence of the relaxation time should be 



T = L 



3 P 1 

l-p2K{L)' 



so that the expression 



p 

T 1 — p 
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2K{L) 



(7) 



(8) 
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FIG. 3: The extrapolated value Huil^oo L'^Ep/(l — p) vs. p. For both boundary dynamics the 
result is consistent with the predicted density dependence. 

should be independent of p. In order to avoid length corrections the limiting case L ^ oo 
is considered by extrapolating the term L'^/r by the algorithm of Bulirsch and Stoer (BST) 
as described in Ref. y 

In the BST algorithm a degree of freedom uj has to be chosen whose size is of the order of 
the correction exponent of the extrapolated series, but whose influence on the result of the 
extrapolation is usually not critical. For the plotted data uj = 1.5 was used, because in this 
region a stable area could be found. Only for the reservoir with p = 0.375 no convergence of 
the algorithm could be achieved. The error of the extrapolation is determined by a variation 
of UJ in the stable area. 

Contrary to expectations the extrapolated values show some dependence on the parameter 
UJ. For UJ = 0.5 which is the assumed correction exponent'' lower values are found but a 
significant p dependence can not be found. 

The results for the different boundary dynamics and particle densities are plotted in Fig. 01 
In consideration of the error bars no density dependence can be found so that we have found 
the correct density dependence. The density dependence found here confirms the dependence 
conjectured by Rubinstein,^^ who proposes a density dependence of r oc z{z — l)N^, where 
2; — 1 is the number of cells which can be entered by an end repton (which is 2d in our 
notation). The factor z — 1 corresponds exactly to p/(l — p) and the additional factor z 
results from the fact that our time scale is 2;-times the time scale used by Rubinstein, so that 
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every time measured in our scale must be multiplied by a factor of z when being compared 
to the Rubinstein scale. We have thus shown that the conjecture is not only valid for the 
Rubinstein model but also for generalizations. 

C. Effective Exponent, details of chain end dynamics 

As in polymer physics the scaling of the viscosity with the polymer weight M is investi- 
gated intensively we will now focus directly on this quantity. Experiments show a scaling of 
ri oc M^'^ for several magnitudes of polymer weight and a crossover to an exponent of 3 in 
the limit of infinitely long chains. 

The effective exponent of the scaling is given by the local slope zat in the log-log-plot, which 
is calculated by 

_ Inrjv+i - lnrjv_i 

ln(iV+l)-ln(iV-l)- ^ ^ 

Fig. m shows the effective exponent for both boundary dynamics with = 3. As abscissa 
values we have chosen 'N'^l'^ because according to Doi the finite size correction should be 
of this order. Compared to the RD model the reservoir-dynamics shifts the effective expo- 
nent to smaller values, with increasing chain length the deviation decreases. While the RD 
model shows an exponent in the region which is found experimentally at least in a small 
range, the values found for the reservoir-dynamics are too small. But using once again 
the BST-algorithm it can be shown that both curves approach zjv = 3 for — > oo: for 
(i = 3 we get for the RD model limjy^oo = 3.03 ± 0.07 and for the reservoir-dynamics 
limjv^oo Z]sf = 3.00 ± 0.03. For other values of d the same qualitative results can be found, 
but for smaller values of d the difference of the curves decreases. 

So we conclude that although the influence of the boundary is much smaller than the in- 
fluence of the parameter d, for some quantities the type of boundary dynamics is crucial. 
We see that the effective exponent is a very sensitive quantity for which the type of model 
which is used for the chain end dynamics is of big influence. 

IV. SUMMARY 

In this paper we have investigated the influence of the boundary dynamics on the viscos- 
ity of polymer melts by comparing the Rubinstein-Duke model to a model with modified 
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FIG. 4: Length dependence of the effective exponent - influence of the boundary dynamics, d = 3 
i.e. p = 6/7. 

boundary dynamics. The model is modified in a way that the entropic force which leads to 
the stretching of the polymer remains the same. In the lattice gas picture this is achieved by 
introducing particle reservoirs at the boundaries which provide the required particle density 
corresponding to the stretching of the polymer. 

We find that the infiuence of particle density and hence the infiuence of the entropic force 
is bigger than the infiuence of the details of the chain end dynamics. The dependence of 
the viscosity on the particle density is quantified in general recovering the conjecture of 
Rubinstein as a special case. Focusing on the scaling of the viscosity we observe that the 
exponent seen in experiments is recovered much better with the Rubinstein-Duke dynamics 
than with its modifications. Nevertheless both models lead to the same exponent of 3 in the 
limit of infinitely long chains. So although the chain end dynamics has smaller infiuence on 
the viscosity of polymer melts than the particle density it nevertheless affects considerably 
the scaling of the viscosity. 
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